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Abstract 


Plasma  bubbles  are  regions  of  depleted  plasma  density  generated  in  the  post¬ 
sunset  equatorial  region  of  the  ionosphere.  Bubbles  significantly  affect  total  electron 
count  (TEC)  and  consequently  alter  communication  and  navigation  capabilities.  Here, 
the  Global  Assimilation  of  Ionospheric  Measurements  Gauss-Markov  (GAIM-GM) 
model  is  studied  in  order  to  assess  its  capability  to  accurately  model  equatorial  plasma 
bubbles.  GAIM-GM  uses  the  Ionospheric  Forecast  Model  (IFM)  as  a  background  state 
modified  through  the  application  of  a  Kalman  Filter  to  incorporate  ionospheric 
observations  such  as  Global  Positioning  System  (GPS)  total  electron  content  (TEC) 
values.  GPS  TEC  values  representative  of  plasma  bubble  conditions  are  modeled  and 
input  into  GAIM-GM  along  with  a  background  IFM  state  and  compared  with  GAIM-GM 
model  runs  with  plasma  bubbles  absent.  The  model  resolution,  bubble  depletion  factor, 
filter  time  constant,  and  geophysical  conditions  are  varied  to  assess  model  response. 
GAIM-GM  is  unable  to  reliably  model  plasma  bubble  magnitude  and  spatial  extent  due 
to  several  characteristics  of  the  model  and  filter.  The  resolution  of  GAIM-GM,  in  both 
regional  and  global  modes,  is  insufficient  to  capture  the  small  scale  of  bubbles. 
Additionally,  the  sun  synchronous  movement  of  perturbations  and  lower  bound  of  the 
Kalman  filter  affect  the  location  and  magnitude  of  depletion  regions  in  the  model. 
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ASSESSMENT  OF  THE  EFFECTS  OF  PLASMA  BUBBLES  ON  GAIM-GM 


I.  Introduction 


1.1  Background 

Plasma  bubbles  are  regions  of  severely  depleted  plasma  densities  commonly 
found  in  the  Earth’s  equatorial  ionosphere.  The  sharp  density  gradients  associated  with 
plasma  bubbles  can  cause  disruption  of  radio  and  satellite  communication.  HF  radio 
communication  and  GPS  are  examples  of  two  systems  the  military  has  grown  dependent 
on  that  are  affected  by  the  ionosphere.  The  presence  of  plasma  bubbles  can  lead  to  radio 
blackouts  and  GPS  signal  outages  that  are  detrimental  to  both  civilian  and  military 
operations. 

The  Global  Assimilation  of  Ionospheric  Measurements  Gauss  Markov  (GAIM- 
GM)  model  is  a  leading  ionospheric  representation  that  is  used  by  the  Air  Force. 
Specifications  and  forecasts  created  by  GAIM-GM  are  used  by  the  Air  Force  in  space  and 
air  operations.  The  accuracy  of  the  model  affects  operational  and  tactical  military 
planning.  Commanders  and  mission  planners  need  to  know  the  reliability  and  tolerances 
of  decision  aids  to  effectively  utilize  resources. 
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1.2  Problem  Statement 


Like  terrestrial  weather  models,  it  is  important  to  understand  the  capabilities  and 
limitations  of  an  ionospheric  model.  Operationally,  it  is  necessary  for  mission  planners 
and  forecasters  to  understand  the  capability  of  GAIM-GM  to  include  plasma  bubbles  in 
its  output.  Computationally,  exploring  the  limits  of  current  systems  forges  the  road  ahead 
for  improvements  to  predictive  capabilities. 

1.3  Research  Objective 

The  research  objective  is  to  determine  the  capability  of  GAIM-GM  to  model 
plasma  bubbles.  GAIM-GM  does  not  include  the  necessary  physics  to  allow  for  the 
consistent  reconstruction  of  plasma  bubble  structures.  Plasma  bubbles  will  be  imposed 
upon  the  input  of  GAIM-GM  and  the  effects  of  the  bubbles  upon  the  output  studied.  The 
results  will  be  useful  to  both  space  weather  forecasters  and  future  model  developers. 
Model  developers  will  be  able  to  incorporate  the  results  into  future  ionospheric  models  to 
better  handle  plasma  bubbles.  Space  weather  forecasters  will  better  understand  the 
capabilities  of  the  present  model  and  be  able  to  tailor  then  products  accordingly. 

1.4  Preview 

The  following  chapters  will  chronicle  the  research  conducted  to  determine  the 
capability  of  GAIM-GM  to  model  plasma  bubbles.  Chapter  2  will  provide  background 
information  on  plasma  bubbles,  GAIM-GM,  and  GPS  ionospheric  measurements.  Next, 
Chapter  3  will  cover  the  methodology  used  to  conduct  the  research,  including  the 
placement  of  plasma  bubble  observations  within  GAIM-GM  and  the  model  parameters 
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adjusted  in  the  study.  Then  results  of  the  tests  involving  various  model  parameters  will 
be  presented  in  Chapter  4.  Finally,  Chapter  5  will  conclude  with  an  assessment  of 
GAIM-GM’s  performance,  recommendations,  and  suggestions  for  future  study. 
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II.  Literature  Review 


This  chapter  will  provide  background  information  on  topics  relevant  to  research 
on  the  capability  of  GAIM-GM  to  model  equatorial  plasma  bubbles.  First,  the  Earth’s 
ionosphere  will  be  discussed  because  plasma  bubbles  occur  in  this  layer  of  the 
atmosphere.  Next,  a  section  describing  plasma  bubbles  will  be  presented,  to  include  the 
formation  process  and  characteristics.  A  discussion  on  calculating  Total  Electron 
Content  (TEC)  from  the  GPS  network  will  follow.  Finally,  GAIM-GM  will  be  described, 
including  discussions  on  the  Ionospheric  Forecast  Model  (IFM)  and  the  Kalman  filter. 

2.1  Ionosphere 

The  ionosphere  is  the  region  of  the  atmosphere  where  significant  numbers  of  free 
thermal  electrons  and  ions  are  present.  Solar  extreme  ultraviolet  (EUV)  radiation  and 
particle  precipitation  are  two  significant  sources  of  energy  input  into  the  ionosphere. 
Within  the  EUV  spectrum  relatively  longer  wavelengths  cause  dissociation  while  shorter 
wavelengths  cause  ionization.  The  sinks  of  ionospheric  energy  are  airglow  and  neutral 
heating  of  the  thermosphere.  The  ion  density  is  controlled  by  the  balance  of 
photoionization  production,  ion  transport,  and  recombination  losses.  The  ionosphere  is 
located  from  60  km  above  the  surface  of  the  Earth  to  altitudes  greater  than  1000  km.  The 
ionosphere  is  generally  divided  into  four  regions,  the  D,  E,  FI,  and  F2  layers,  as  shown  in 
Figure  1.  The  F2  layer  has  the  greatest  ion  density  and  consequently  impact  on  radio  and 
satellite  signal  propagation.  At  equatorial  latitudes  the  plasma  density  of  the  ionosphere 
is  higher  during  the  day  due  to  solar  photoionization,  while  at  night  solar  photoionization 
ends  and  recombination  significantly  decreases  the  plasma  densities  of  the  D  and  E 
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0) 


Figure  1.  The  layers  of  the  Earth's  daytime  ionosphere.  From  Schunk  and  Nagy,  2009. 

layers.  Plasma  densities  in  the  ionosphere  vary  with  season,  time  of  day,  solar  cycle,  and 
geomagnetic  activity.  Irregularities  in  the  plasma  density  are  known  to  occur  and  some 
of  the  most  severe  ones  are  called  plasma  bubbles  (Schunk  and  Nagy,  2009). 

2.2  Plasma  Bubbles 

Plasma  bubbles  are  regions  of  reduced  F  region  plasma  densities  that  form  along 
the  bottomside  F  region  (200  km)  and  extend  vertically  into  the  topside  ionosphere.  Also 
known  as  plasma  depletions  or  holes,  plasma  bubbles  have  characteristic  spatial  and 
temporal  patterns  that  have  been  studied  for  more  than  30  years  (Sinha  and  Raizada, 
2000).  The  bubbles  increase  ionospheric  scintillations  resulting  in  disruption  or 
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degradation  of  transionospheric  communications  and  navigation  systems  (Sahai  et  al, 

2000). 

Plasma  bubbles  were  first  observed  as  “ion  density  biteouts”  by  Hanson  and 
Sanatai  in  1973  using  in  situ  satellite  probes.  The  bubbles  manifest  themselves  as  regions 
of  low  intensity  in  night  airglow  images.  Modern  methods  of  detecting  these  ionospheric 
anomalies  include  ionosondes,  incoherent  scatter  radars,  in  situ  rockets  and  satellites,  and 
wide  angle  airglow  imaging  systems  (Sinha  and  Raizada,  2000).  Currently,  GPS  offers 
continuous  monitoring  of  the  generation  and  evolution  of  plasma  bubbles  on  a  large 
scale,  due  to  its  permanent  operation  and  near  global  coverage  (Portillo  et  al,  2008). 

2.2.1  Development 

Plasma  bubbles  begin  as  Rayleigh-Taylor  instabilities  in  the  bottomside  of  the 
equatorial  F  layer  (Burke  et  al,  2004).  The  bubbles  form  as  density  perturbations  which 
are  believed  to  be  caused  by  gravity  waves.  The  resulting  instability  propagates  upwards 
through  the  post  sunset  F  region  (Li  et  al,  2009). 

Plasma  bubbles  are  the  result  of  dynamo  electric  fields,  which  are  generally 
eastward  during  the  day  and  westward  at  night.  During  the  day  the  eastward  electric  field 
produces  an  upward  ExB  drift  of  the  F  region  plasma.  Near  dusk  the  eastward  electric 
field  is  enhanced,  referred  to  as  the  prereversal  enhancement,  before  it  turns  westward. 
The  F  region  consequently  rises  due  to  the  stronger  upward  ExB  drift  as  the  ionosphere 
corotates  into  the  darkness.  In  the  darkness  the  lower  part  of  the  F  region  is  further 
depleted  by  recombination  where  the  ion-neutral  collision  frequency  is  high.  The  result 
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is  a  raised  F  region  with  a  steep  vertical,  bottomside  density  gradient  (Schunk  and 
Demars,  2003). 

The  steep  density  gradient,  shown  at  the  top  of  Figure  2,  is  supported  by  the 
strong  magnetic  field.  Initially  no  electric  fields  are  present  but  when  a  perturbation 
occurs,  charge  build  up  leads  to  electric  fields  in  the  y  component.  Electric  fields  in  the 
positive  y  direction  lead  to  an  upwards  ExB  drift  that  allows  the  instability  to  grow  into 
small  density  irregularities,  or  spread  F.  During  severe  spread  F,  plasma  bubbles  form  on 
the  bottomside  F  layer  and  drift  upwards  with  speeds  that  generally  vary  from  100  to  500 
m/s  (Schunk  and  Nagy,  2009). 

As  the  instability  moves  upwards  the  entire  north-south  extent  of  the  magnetic 
flux  tube  becomes  depleted.  The  bubble  takes  the  form  of  a  vertically  elongated  wedge 
of  depleted  plasma,  resembling  an  upside  down  watermelon  slice,  that  can  be  hundreds  of 
kilometers  wide  (Schunk  and  Demars,  2003). 

2.2.2  Characteristics 

Plasma  bubbles  are  most  likely  to  form  in  regions  where  the  ends  of  magnetic  flux 
tubes  encounter  the  dusk  terminator  simultaneously.  At  longitudes  where  the  magnetic 
declination  is  zero,  plasma  bubble  formation  is  favored  near  the  equinoxes.  In  regions  of 
westward  (eastward)  magnetic  declination,  the  minimum  angle  between  the  flux  tubes 
and  the  dusk  terminator  shifts  towards  the  December  (June)  solstice.  Most  plasma 
bubbles  are  observed  within  15  degrees  of  the  geomagnetic  equator  (Burke  et  al,  2004). 
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Equilibrium  state 
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Perturbed  state 


Figure  2.  Rayleigh-Taylor  Instability  configuration.  X  points  to  the  south,  y  points  to  the  east,  and  z 
points  upward.  Adapted  from  Schunk  and  Nagy,  2009. 

Plasma  bubbles  typically  form  at  2000  LT  and  remain  as  a  density  depleted  region 
through  2400  LT.  Bubbles  can  also  last  throughout  the  night  until  0600  LT  when 
sunlight  replenishes  the  plasma  in  the  F  layer  (Li  et  al,  2009).  Plasma  within  the  bubble 
is  typically  depleted  by  a  factor  of  about  10,  but  factors  of  100  and  1000  and  also  fairly 
common  (Schunk  and  Demars,  2003). 

Plasma  bubbles  have  been  observed  to  reach  altitudes  greater  than  1500  km  and 
have  zonal  widths  of  several  hundred  kilometers  (Sahai  et  al,  2000).  The  zonal  distance 
between  bubbles  is  on  the  order  of  the  bubble  width  and  has  been  recorded  between  400 
and  1600  km.  Plasma  bubbles  have  been  shown  to  have  both  eastward  and  westward 
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tilts,  up  to  15  degrees,  with  height  (Sinha  and  Raizada  2000).  Additionally,  bubbles  have 
been  recorded  corotating  with  the  Earth  while  others  super-rotate  with  speeds  up  to  3 
km/s  (Schunk  and  Demars,  2003). 

Table  1  summarizes  plasma  bubble  characteristics. 

2.3  Inferring  TEC  from  GPS  Measurements 

Data  collected  from  a  continuously  operating  global  network  of  Global 
Positioning  System  (GPS)  receivers  can  be  used  to  produce  global  maps  of  ionospheric 
total  electron  content.  The  GPS  network  consists  of  hundreds  of  dual  frequency  geodetic 
quality  receivers  and  31  GPS  satellites.  Each  GPS  receiver  tracks  between  4  and  12 
satellites  simultaneously  over  a  site.  The  TEC  in  the  line  of  site  between  the  ground 
station  and  satellite  can  be  computed  using  information  carried  in  the  GPS  signal 
(Mannucci  et  al,  1998). 

The  GPS  receivers  in  the  global  network  generate  two  delay  observables  for  each 
satellite  being  tracked,  pseudorange  delay  and  carrier  phase  delay.  Each  observable  is 
acquired  at  the  two  GPS  transmission  frequencies,  fl  (1.5754  GHz)  and  f2  (1.2276  GHz). 


Table  1.  Plasma  Bubble  Characteristics.  Adapted  from  Nava,  2011. 


Characteristic 

Typical  Values 

Depletion  Extent 

Factor  of  10,  100,  1000 

N-S  Extent 

Several  thousands  of  km 

Width 

Several  hundreds  of  km 

Apex  Height 

>  1500  km 

Start  Time 

2000  LT 

Dissipation  Time 

0600  LT  (sunrise) 

Upward  Drift 

100  to  500  m/s  (slow) 

500  m/s  to  5  km/s  (fast) 

Eastward  Drift 

40-190  m/s 

Tilt 

10-15  degrees  to  the  west 

9 


The  frequency  differenced  phase  delays  provide  very  precise  measurements  of  TEC 
changes  but  contain  an  overall  bias  associated  with  integer  cycle  ambiguities.  The 
frequency  difference  pseudorange  data  provide  an  absolute  measure  of  the  total 
ionospheric  delay  between  the  satellite  and  receiver  but  have  significantly  more  noise 
than  the  phase  based  data.  In  addition,  biases  are  created  by  delays  in  the  satellite  and 
receiver  hardware  (Mannucci  et  al,  1998). 

A  TEC  observable  is  formed  from  the  frequency  differenced  phase  delay  by 
adjusting  its  level  based  on  the  pseudorange  delay.  The  TEC  observable  for  each 
measurement  is  the  phase  based  ionospheric  delay  added  to  the  bias.  One  nanosecond  of 
delay  corresponds  to  2.85  TEC  units  (Mannucci  et  al,  1998). 

Vertical  TEC  is  produced  from  a  slant  TEC  observation  measured  from  the 
satellite  to  the  ground  station.  The  calculation  utilizes  a  fictitious  ionospheric  pierce 
point  defined  by  the  intersection  of  the  satellite-receiver  line  of  sight  with  a  thin  spherical 
shell  surrounding  the  Earth,  often  at  an  altitude  of  450  km.  The  location  of  this  point  is 
then  dropped  vertically  in  altitude  to  the  Earth’s  surface  where  it  is  assigned  a  latitude- 
longitude  location.  The  latitude-longitude  location  of  the  vertical  TEC  observation  will 
necessarily  differ  from  the  location  of  the  ground  station  for  all  elevation  angles  less  than 
90  degrees.  Figure  3  shows  the  maximum  radial  distance  possible  for  a  vertical  TEC 
measurement  from  the  associated  ground  station  for  several  locations  worldwide 
(Mannucci  et  al,  1998). 
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Figure  3.  A  distribution  of  GPS  ground  stations.  Ellipses  surrounding  each  ground  station  represent 
the  maximum  distance  a  vertical  TEC  observation  will  be  placed  from  the  ground  station,  assuming 
an  ionospheric  shell  pierce  point  of  450  km.  From  Howe  and  Runciman,  1998. 


2.4  Global  Assimilation  of  Ionospheric  Measurements  (GAIM)  Model 

GAIM  is  a  physics-based  data  assimilation  model  of  the  ionosphere  developed  by 
the  Center  for  Atmospheric  and  Space  Sciences  at  Utah  State  University  (USU).  At 
present  the  USU  team  has  fielded  GAIM  Gauss-Markov  (GAIM-GM),  using  the 
Ionospheric  Forecast  Mode  (IFM)  as  a  background  model,  as  the  operational  Air  Force 
ionospheric  model  at  the  Air  Force  Weather  Agency.  USU  is  currently  developing 
GAIM  Full  Physics  (GAIM-FP)  which  uses  the  Ionosphere-Plasmasphere  Model. 
GAIM-GM  utilizes  a  Gauss-Markov  Kalman  filter  to  ingest  ionospheric  observations  and 
weight  the  observations  against  the  IFM  background  density.  Ionospheric  observations 
from  several  different  sources  are  able  to  be  assimilated  and  GAIM-GM  can  be  expanded 
to  include  a  larger  network  of  measurements  (Schunk  et  al,  2004). 
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GAIM-GM  combines  a  physics  based  model,  the  IFM,  with  ionospheric 
observations.  The  ionospheric  observations  can  include  in  situ  electron  densities  from 
satellites,  bottomside  electron  density  profiles,  GPS  TEC  measurements,  satellite 
occultations,  and  nighttime  UV  radiances.  Of  these  observations,  GPS  TEC 
measurements  and  ionosondes  are  the  most  commonly  used,  with  the  GPS  TEC 
measurements  representing  the  bulk  of  the  ingested  data.  A  Kalman  filter  is  used  to 
assign  weights  to  the  relative  accuracies  of  the  observations  and  model  inputs.  The  data  is 
then  processed  to  reconstruct  an  output  of  the  ionosphere.  The  outputs  of  GAIM-GM 
include  a  three  dimensional  electron  density,  TEC,  NmF2,  and  HmF2.  GAIM-GM 
produces  these  outputs  in  user  specified  global  or  regional  windows  with  a  grid  system  of 
44  nodes  meridionally  and  24  nodes  zonally.  In  the  global  mode  resolution  is  4.67 
degrees  of  latitude  by  15  degrees  of  longitude  while  the  finest  regional  resolution  is  1 
degree  of  latitude  by  3.75  degrees  of  longitude.  Boundary  conditions  for  the  regional  run 
of  GAIM-GM  are  set  from  the  output  of  a  prior  global  GAIM-GM  model  run  with  the 
same  geophysical  conditions  (Schunk  et  al,  2004). 

GAIM-GM  uses  calculated  GPS  satellite  and  ground  station  biases  to  correct  GPS 
TEC  measurements.  Each  GPS  satellite  and  ground  station  has  a  known  bias  measured  in 
nanoseconds,  resulting  from  satellite  drag  and  clock  differences.  Also,  an  error  is 
calculated  and  associated  with  each  satellite  and  ground  station  regarding  the  precision  of 
the  measurement  (Schunk  et  al,  2004). 

Additionally  GAIM-GM  includes  several  parameters  that  can  be  set  by  the  user 
prior  to  each  run.  The  time  constant  is  the  relaxation  time  the  model  uses  to  decay 
perturbations  to  zero.  The  time  constant  can  be  changed  along  an  integer  range  from  1  to 
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12  hours  and  the  default  is  5  hours.  GAIM-GM  moves  perturbations  sun  synchronously 
and  the  effect  of  lowering  the  time  constant  is  to  reduce  to  westward  extent  the 
perturbation  field  will  take  to  relax  to  zero.  The  latitude  correlation  describes  the 
meridional  distance  the  perturbation  field  is  allowed  to  extend.  The  latitude  correlation 
can  be  adjusted  along  integer  values  from  2  to  10  degrees  and  the  default  value  is  10 
degrees.  Decreasing  the  latitude  correlation  lowers  north-south  distance  of  the 
perturbation  field  (Scherliess  et  al,  2006). 

2.4.1  IFM 

The  IFM  is  a  physics-based  numerical  model  of  the  global  ionosphere  that 
calculates  three  dimensional,  time  dependent  density  distributions  for  four  major  ions 
(NO+,  O2+,  N2+,  0+)  at  E  region  altitudes,  two  major  (O+,  NO+)  and  two  minor  (N2+, 
02+)  ions  at  F  region  altitudes.  Additionally  the  IFM  calculates  ion  and  electron 
temperatures  at  both  E  and  F  region  altitudes  as  well  as  H+  densities  in  the  F  region  and 
topside  ionosphere.  The  primary  output  of  the  IFM  is  a  three  dimensional  electron 
density  ranging  in  altitude  from  90  to  1600  km  with  a  resolution  of  4  km  in  the  E  region 
and  20  km  in  the  F  region.  The  output  density  has  a  3  degree  latitude  resolution  and  a 
7.5  degree  longitude  resolution  (Schunk  et  al,  2004). 

IFM  inputs  include  FI 0.7,  year,  day,  start  time,  duration  of  model  run,  and  the 
temporal  variation  of  Kp  from  3  hours  prior  to  the  start  time  through  the  end  of  the  model 
run.  The  IFM  is  based  on  a  numerical  solution  of  the  ion  and  electron  continuity, 
momentum,  and  energy  equations.  The  model  is  able  to  account  for  the  following 
physical  processes:  1)  field-aligned  diffusion  due  to  density  and  temperature  gradients, 
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gravity,  and  the  ambipolar  electric  field;  2)  cross-field  electrodynamic  drifts  due  to  both 
magneto  spheric  and  dynamo  electric  fields;  3)  ion  production  due  to  UV  and  EUV  solar 
radiation,  starlight,  and  auroral  electron  precipitation;  4)  numerous  energy-dependent 
chemical  reactions;  5)  thermospheric  winds;  6)  neutral  composition  changes;  7)  a  myriad 
of  elastic  and  inelastic  heating  and  cooling  processes.  Additionally  IFM  accounts  for  the 
offset  between  the  geographic  and  geomagnetic  poles  (Schunk  et  al,  2004). 

2.4.2  Kalman  Filter 

Howe  and  Runcimen  demonstrated  the  practicality  of  using  GPS  TEC 
measurements  within  the  context  of  a  Kalman  filter  to  produce  four  dimensional  maps  of 
the  ionosphere  (1998).  The  Kalman  filter  is  a  well  established  least  squares  procedure. 
The  filter  provides  an  optimal  method  for  assimilating  many  different  data  types  into  a 
time-dependent  model.  The  specification  of  error  covariances  for  the  model  and 
observation  error  is  critically  important  to  the  proper  functioning  of  the  filter.  The 
Kalman  filter  utilizes  a  series  of  six  equations  to  minimize  the  error  at  each  time  step.  A 
forecast  is  made  of  the  model  state,  the  ionosphere,  using  Equation  1.  The  state  vector  x 
is  made  by  linearly  relating  the  state  vector  at  a  time  t  +  1  to  the  state  at  the  previous  time 
step  t  through  the  transition  matrix  M,  where  q  is  in  the  model  error.  The  model  state  at 
the  new  time  step  is  linearly  related  to  the  model  state  at  the  previous  time  step  through  a 
transition  matrix  and  a  forcing  term.  A  forecast  model  state  estimate  comes  with  the 
forecast  model  state  error  covariance,  Equation  2,  which  depends  on  the  previous  error 
and  the  model  error.  All  information  about  the  model  state  error  is  carried  forward  from 
one  time  step  to  the  next  by  the  model  error  covariance  matrix  P  and  no  information  older 
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than  one  time  step  is  explicitly  retained.  The  transition  model  error  covariance  is  given 
by  matrices  Q  and  R.  Measurements  y°  are  represented  by  Equation  3,  where  H  is  the 
measurement  matrix  and  s  is  the  observation  error.  The  model  state  forecast  is  combined 
with  the  data  using  the  Kalman  gain  matrix  K,  given  in  Equation  4,  to  produce  the 
analysis  model  estimate,  Equation  5,  which  is  the  optimal  model  state  estimate  given  a 
knowledge  of  the  model  state  and  the  data  and  model  error  covariances.  Equation  6 
produces  the  analysis  model  error  covariance  which  is  reduced  from  the  forecast  error  by 
an  amount  related  to  the  resolving  power  of  the  data  (Howe  and  Runcimen,  1998). 


=  Lx  +  T] 

(1) 

pf  =  LPLt  +  Q 

(2) 

y°  =  Hx  +  £ 

(3) 

pfHT(HPfHT  +  R)1 

(4) 

=  x?  +  K(y°  —  Hx? ) 

(5) 

Pa  =  (/  -  KH)pf 

(6) 

GAIM-GM  uses  a  Gauss-Markov  Kalman  filter  to  assimilate  ionospheric 
observations  with  IFM  background  states.  The  background  ionospheric  density  is 
obtained  from  the  IFM’s  electron  density  output  while  ionospheric  measurements  are 
assigned  to  the  perturbation  density.  The  Kalman  filter  applies  the  calculated 
perturbation  density  to  the  background  density  to  output  a  solution  for  each  grid  location, 
according  to  Equations  7  and  8. 

N  —  Nifm  +  Npert  (7) 
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(8) 


Nt+1  —  I  Nt 
‘'pert  ‘-'‘'pert 


A  transition  matrix  L  is  used  by  the  filter  at  each  time  step  to  propagate  perturbations 
forward  through  time,  described  in  Equation  8.  The  transition  matrix  has  translation  and 
decay  components  included.  Perturbations  move  forward  in  a  magnetic  sun-synchronous 
frame  and  decay  towards  a  value  of  zero  according  to  a  relaxation  parameter  (Scherliess 
Et  AL,  2006). 


The  Kalman  filter  also  incorporates  a  model  error  covariance  matrix,  used  to 
estimate  the  uncertainties  in  the  model  state  which  is  the  perturbation  field.  The  model 
error  covariance  matrix  has  been  predetermined  from  1107  individual  two-day  runs  of  the 
IFM  that  spanned  a  range  of  ionospheric  forcing  conditions.  The  filter  also  applies 
correlations  in  both  the  zonal  and  meridional  directions.  Zonal  correlation  are  tied  to  the 
time  constant  while  meridional  correlations  are  set  by  the  latitude  correlation  (Scherliess 
Et  AL,  2006). 


The  five  types  of  observations  the  Kalman  filter  in  GAIM-GM  can  ingest  are 
TECs  from  GPS,  ionosondes,  in  situ  electron  densities,  electron  density  profiles  from 
ground  based  radars,  and  optical  emission  profiles.  The  most  common  data  types  used  by 
GAIM-GM  are  GPS  TEC  measurements  and  ionosondes  (Thompson  et  al,  2006). 
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III.  Methodology 


This  chapter  describes  the  methodology  used  in  conducting  the  research.  The  first 
section  discusses  the  plasma  bubble  structure  and  placement.  Then  the  calculation  of 
slant  TEC  and  generation  of  GPS  files  is  covered.  The  next  section  details  adjustments 
and  considerations  within  GAIM-GM  to  include  station  and  satellite  biases,  the 
plasmasphere,  the  time  constant,  and  the  latitude  correlation.  Finally,  the  conditions 
GAIM-GM  was  tested  under  are  outlined. 

3.1  Generating  Plasma  Bubble  Observations 

First,  a  method  to  input  plasma  bubble  observations  into  GAIM-GM  was 
developed.  GAIM-GM  uses  the  IFM  background  state  modified  by  GPS  slant  TEC  data. 
The  IFM  is  a  physics-based  model  that  accounts  for  the  date,  solar,  and  geomagnetic 
conditions.  The  IFM  does  not  use  ionospheric  observations  and  does  not  contain  the 
necessary  physics  to  model  plasma  bubbles.  GPS  slant  TEC  data  are  measurements  of 
the  ionosphere  calculated  from  the  phase  delay  of  signals  from  GPS  satellites.  Of  these 
two  data  sources  into  GAIM-GM,  real-world  plasma  bubble  observations  would  be 
present  in  the  GPS  slant  TEC  data.  Consequently  utilities  were  developed  to  place 
plasma  bubble  observations  into  synthetic  GPS  slant  TEC  data. 

In  order  to  isolate  the  effect  of  plasma  bubbles  on  GAIM-GM  output,  it  is 
necessary  to  have  consistency  between  the  IFM  and  GPS  slant  TEC  inputs.  If  the  GPS 
measurements  differ  from  the  IFM  background  state  the  Kalman  filter  will  introduce 
perturbation  fields  to  reconcile  the  differences.  The  perturbation  fields  have  far  reaching 
effects  in  the  model  grid  due  to  the  time  constant  and  latitude  correlation,  hindering  any 
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plasma  bubble  signature.  Therefore  the  IFM  electron  density  field  was  used  as  the 
background  from  which  to  construct  artificial  GPS  slant  TEC  data.  Slant  TECs  were 
calculated  through  the  IFM  electron  density  both  with  and  without  bubbles  present.  The 
slant  TECs  calculated  with  bubbles  absent  created  a  set  of  measurements  identical  to  the 
IFM  background  state  used  by  GAIM-GM,  limiting  the  presence  of  perturbation  fields. 
GAIM-GM  was  then  run  using  slant  TEC  measurements  with  bubbles  present  to  be 
evaluated  against  the  bubbles  absent  case.  Both  IFM  and  slant  TEC  files  were  provided 
to  GAIM-GM  in  15  minute  increments  for  each  model  day. 

3.1.1  Plasma  Bubble  Structure 

Plasma  bubbles  were  imposed  upon  the  IFM  electron  density  field  according  the 
method  described  by  Nava  (2011).  First,  the  original  resolution  of  the  IFM  electron 
density  grid  was  doubled  to  produce  a  grid  with  3.75  degree  longitudinal  and  1.5  degree 
latitudinal  spacing.  The  variable  spacing  between  the  altitude  grid  points  was  halved  in 
order  to  double  the  vertical  resolution.  The  electron  densities  were  linearly  interpolated 
to  determine  their  values  at  the  newly  created  IFM  grid  points.  The  grid  resolution  was 
increased  to  create  plasma  bubbles  of  a  realistic  size.  Figure  4  prescribes  the  method  to 
measure  bubble  width  in  the  modified  IFM  electron  density  grid.  Because  slant  TECs 
calculated  through  the  grid  take  an  average  of  the  nearest  eight  grid  points  in  three 
dimensional  space,  the  effective  width  of  the  bubble  spans  from  one  undepleted  column 
of  grid  points,  across  the  depleted  column,  and  to  the  next  undepleted  column  of  grid 
points.  For  the  modified  IFM  grid  with  3.75  degrees  between  grid  points  the  result  is  a 
bubble  width  of  7.5  degrees,  a  large  width  for  a  plasma  bubble  but  still  physically 
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Figure  4.  Depletion  width  of  a  plasma  bubble  in  the  modified  IFM  electron  density  resolution. 
Adapted  from  Nava,  2011. 

reasonable.  In  the  original  7.5  degree  longitudinal  resolution  of  the  IFM,  a  bubble  width 
would  be  15  degrees,  too  wide  to  be  physically  reasonable  and  necessitating  the 
modification  to  increase  the  IFM  resolution  (Nava  2011). 

Next,  specific  grid  points  were  chosen  to  form  the  outline  of  the  plasma  bubble. 
The  bubble  was  shaped  like  an  upside  down  watermelon  slice  elongated  along  Earth’s 
magnetic  field  lines.  The  base  meridional  length  of  the  bubble  was  specified,  from  which 
the  apex  height  of  the  bubble  is  calculated  from  Equation  3, 

r  =  R0sin26  (3) 

where  r  is  the  dipole  field  line  radius,  9  is  the  colatitude,  and  Ro  is  the  radial  distance  to 
the  point  where  the  dipole  field  crosses  the  geomagnetic  equator.  The  electron  density  at 
each  chosen  grid  point  was  then  divided  by  the  specified  depletion  factor.  A  side  view  of 
a  plasma  bubble  in  the  IFM  electron  density  is  shown  in  Figure  5  (Nava  2011). 
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Figure  5.  Longitudinal  view  of  a  plasma  bubble  in  the  IFM  electron  density  grid.  Scale  is  in  units  of 
electrons  per  square  meter. 

3.1.2  Plasma  Bubble  Placement 

Although  plasma  bubbles  have  been  observed  globally  across  the  equatorial 
regions  of  the  Earth  it  was  decided  to  focus  on  bubbles  occurring  over  the  South 
American  sector.  South  America  contains  the  highest  density  of  equatorial  ground 
stations  and  provides  a  distribution  that  has  a  large  latitudinal  extent.  Plasma  bubbles 
were  imposed  on  the  IFM  electron  density  grid  over  South  America.  The  bubbles  were 
formed  at  2000  LT  and  removed  at  0600  LT  to  match  typical  characteristics.  Because  the 
bubble  upward  velocity  during  its  formation  can  be  large,  the  15  minute  time  resolution 
of  the  IFM  is  less  than  the  time  it  takes  a  bubble  to  rise  from  200  to  1400  km,  therefore  a 
bubble  appears  in  its  full  form  from  one  IFM  time  step  to  the  next.  Bubbles  were  placed 
7.5  degrees  apart  as  measured  from  edge  to  edge  and  shown  in  Figure  4.  When  a 


20 


sequence  of  bubbles  was  used,  each  bubble  would  form  at  2000  LT  for  its  longitude  and 
the  next  bubble  would  form  7.5  degrees  to  the  west  one  hour  later.  Once  a  bubble  has 
formed,  three  different  methods  of  movement  were  modeled  to  test  GAIM-GM.  The  first 
method  kept  the  bubbles  corotating  with  the  Earth,  a  common  characteristic  of  bubbles, 
shown  in  Figure  6.  Some  bubbles  super-rotate  which  led  to  a  second  condition  where  the 
bubbles  would  form  at  2000  LT  and  then  allowed  to  super-rotate  at  230  m/s,  or  7.5 
degrees/hr,  depicted  in  Figure  7.  Figure  8  is  the  third  and  unrealistic  case,  with  the 
bubble  moving  sun  synchronously  to  highlight  the  performance  of  the  time  constant. 

3.1.3  Slant  TEC  calculation 

With  the  depleted  region  representing  the  plasma  bubble  placed  in  the  electron 
density  grid  the  following  step  was  performed  to  calculate  slant  TECs  through  the 
artificial  ionosphere.  The  method  for  calculating  slant  TEC  was  described  by  Nava  and 
is  an  adaptation  of  the  code  generated  by  USU.  The  code  creates  a  linear  path  between 
the  satellite  and  ground  station  and  places  this  path  in  the  context  of  the  IFM  electron 
density  grid.  Points  are  selected  every  5  km  along  this  path  where  the  8  nearest  grid 
points,  three  dimensionally,  are  linearly  interpolated  to  construct  a  representative  electron 
density  for  the  path  point.  The  path  points  are  then  summed  and  a  conversion  factor  is 
applied  to  change  from  electron  density  to  the  path  integrated  TEC.  The  calculated  TEC 
thus  represents  the  TEC  through  the  IFM  electron  density  grid  along  the  path  from  the 
satellite  to  the  ground  station  (Nava  2011). 
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Figure  6.  Formation  of  4  corotating  plasma  bubbles  for  day  82,  F10.7  =  113  SFUs.  Bubbles  are 
placed  in  modified  IFM  electron  density  grid  and  numbered  according  to  formation  order.  Units  are 
in  TEC  electrons  per  square  meter. 
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Figure  7.  Evolution  of  4  super-rotating  bubbles  for  day  82,  F10.7  =  113  SFUs.  Bubbles  placed  in 
modified  IFM  electron  density  grid  and  numbered  according  to  order  while  moving  east  at  230  m/s. 
Units  are  in  TEC  electrons  per  square  meter. 
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Figure  8.  Evolution  of  1  sun  synchronous  bubble  for  day  82,  F10.7  =  113  SFUs.  Bubble  placed  in 
modified  IFM  electron  density  grid.  Units  are  in  TEC  electrons  per  square  meter. 
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3.1.4  GPS  Input  Files 

GPS  files  are  ingested  by  GAIM-GM  and  reduced  to  a  readable  format  for  the 
model.  Each  GPS  file  contains  the  data  from  one  ground  station  for  a  one  hour  time 
period.  The  data  is  organized  into  22  columns  which  include  time,  date,  satellite 
identification  information,  look  angle  data,  derived  TEC  data,  satellite  location,  and 
ground  station  location.  GAIM-GM  uses  the  date,  time,  satellite  identification,  satellite 
location,  ground  station  location,  TEC,  and  TEC  error  when  ingesting  GPS  observations. 
New  GPS  files  were  written  as  text  files  in  MATLAB  into  the  format  expected  by 
GAIM-GM.  Time,  date,  satellite  identification,  satellite  location,  ground  station  location, 
TEC,  and  TEC  error  were  all  specified  for  each  GPS  file.  The  TEC  error  was  set  at  3 
TECU  to  allow  for  differences  in  representativeness  of  a  given  slant  path  through  the 
modified  electron  density  grid.  Finally,  the  computed  TEC  value  from  the  modified 
electron  density  grid  was  written  to  the  file.  GPS  input  files  were  created  every  hour  for 
each  ground  station.  Each  ground  station  observes  between  5  and  10  satellites  at  each 
time  step  with  position  and  TEC  calculated  at  15  minute  increments.  GPS  files  were 
translated  into  reduced  .gps  format  that  is  readable  by  GAIM-GM.  The  files  were  then 
transferred  to  the  GAIM-GM  GPS  directory  to  be  ingested  into  the  model. 

3.2  GAIM  Adjustments 

3.2.1  Satellite  and  Station  Biases 

GAIM-GM  corrects  the  GPS  TEC  measurements  for  known  satellite  and  ground 
station  biases.  Under  normal  operation  the  satellite  and  ground  station  biases  are 
downloaded  from  an  archive  at  the  University  of  Bern  in  Switzerland.  Each  file  contains 


25 


the  monthly  bias  and  root  mean  square  uncertainty  of  the  bias,  both  expressed  in 
nanoseconds  (GAIM  2010).  Because  this  research  is  using  synthetic  GPS  files,  historic 
satellite  and  ground  station  biases  do  not  apply.  In  real  GPS  measurements  it  is  necessary 
to  have  a  satellite  and  ground  station  bias  adjustment  to  account  for  measurement  error, 
but  since  the  calculation  through  the  IFM’s  electron  density  does  not  use  equipment  to 
observe  the  TEC,  the  measurement  is  exact  and  the  satellite  and  ground  station  biases  are 
zero.  The  slant  TECs  calculated  are  an  exact  solution  because  the  IFM  creation  of  the 
ionosphere  used  to  derive  the  TEC  measurement  is  a  known  quantity.  For  the  purposes 
of  this  study  both  the  satellite  and  ground  station  biases  were  set  at  a  nominal  value  of 
0.001  ns,  the  lowest  value  accepted  by  the  model. 

3.2.2  Plasmasphere 

In  addition  to  the  satellite  and  ground  station  biases,  GAIM-GM  accounts  for  the 
TEC  contribution  from  the  plasmasphere  in  its  use  of  GPS  files.  Because  GAIM-GM 
calculates  an  electron  density  from  90  to  1380  km  and  satellite  observations  measure  the 
entire  TEC  from  satellite  to  ground  station,  TEC  from  1380  km  and  upwards  in  the 
plasmasphere  must  be  subtracted  from  the  GPS  TEC  observation  to  provide  GAIM-GM 
with  a  TEC  specific  to  its  altitudes.  Again,  in  this  study  the  slant  TECs  used  in  the  GPS 
files  were  generated  from  the  modified  IFM  electron  density  grid  from  90  to  1598  km. 
Thus  there  is  no  plasmasphere  contribution  in  the  GPS  file  TEC  measurement  and  no 
need  to  subtract  this  quantity  from  the  TEC  contained  in  the  GPS  input  files.  USU 
provided  a  flag  in  the  GAIM-GM  startup  file  that  removed  the  subtraction  of  the 
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plasmasphere  within  GAIM-GM,  which  allows  the  Kalman  filter  to  receive  the 
appropriate  TEC  value  for  this  study. 

The  slant  TEC  path  through  the  IFM  is  calculated  for  IFM  altitudes  from  90  to 
1598  km.  The  top  altitude  of  GAIM  is  1380  km,  leaving  a  218  km  difference  between 
the  tops  of  the  two  models.  When  GAIM-GM  makes  its  usual  plasmasphere  correction 
the  model  subtracts  the  electron  density  contribution  from  1380  km  up  to  the  satellite 
orbit.  Slant  TECs  used  in  this  study  were  derived  from  the  IFM,  meaning  an  extra 
contribution  of  electron  density  was  added  from  1380  km  to  1598  km.  As  part  of  this 
research,  a  study  quantified  the  TEC  contribution  from  1380  to  1598  km  in  the  IFM.  It 
was  found  that  this  218  km  layer  contributes  a  maximum  value  of  0.23  TECU  in  the 
dayside,  midlatitude  ionosphere  and  0.05  TECU  in  the  nightside,  equatorial  ionosphere. 
Because  the  TEC  contribution  from  this  layer  is  small  compared  to  the  background  TEC 
values,  the  effect  is  assumed  to  be  negligible  and  was  not  corrected. 

3.3  Running  the  Models 

The  IFM  is  an  input  to  GAIM-GM  and  must  be  run  prior  to  GAIM-GM.  The  IFM 
requires  a  24  hour  model  time  warm  up  prior  to  producing  useful  output.  The  IFM  was 
run  for  two  day  series’  to  construct  a  useable  input  for  GAIM-GM.  The  first  day  of  each 
series  was  discarded  as  a  warm  up  day  and  the  second  day  was  used  for  input  to  GAIM- 
GM.  For  example,  to  create  useable  IFM  output  for  GAIM-GM  for  days  8 1  and  82,  2- 
two  day  runs  of  the  IFM  were  completed.  The  first  IFM  run  was  for  days  80  and  81  and 
the  second  IFM  run  was  for  days  81  and  82.  In  each  of  these  runs  the  second  day  of 
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output  was  used  as  the  input  to  GAIM-GM  making  a  useable  set  of  IFM  files  for  days  8 1 
and  82. 

GAIM-GM  also  requires  a  warm  up  day  to  create  correct  ionospheric  output.  For 
example,  to  create  output  for  day  82  GAIM-GM  was  run  for  days  81  and  82.  This 
requires  IFM  and  GPS  input  for  both  days.  For  runs  containing  bubble  observations  in 
the  GPS  files,  the  second  day  of  input  files  contained  slant  TEC  values  from  the  IFM 
with  bubbles  present  while  the  first  day  was  composed  of  slant  TEC  values  from  the  IFM 
without  bubbles.  For  cases  involving  a  regional  GAIM-GM  the  model  must  previously 
be  run  globally  to  create  boundary  conditions.  The  global  output  boundary  conditions 
were  created  from  GAIM-GM  without  GPS  inputs  present. 

3.4  Results  Analysis  Format 

The  results  in  this  paper  are  presented  three  formats:  GAIM-GM  TEC  output, 
percent  difference  maps,  and  skill  scores.  The  GAIM-GM  TEC  output  is  a  display  of  the 
TEC  in  longitude  and  latitude  coordinates.  The  percent  difference  maps  show  a 
quantitative  representation  of  the  spatial  extent  of  the  plasma  bubble  depletions  in  GAIM- 
GM.  The  skill  scores  yield  a  numerical  representation  of  the  relative  accuracy  of  the 
GAIM-GM  representation  compared  to  climatology. 

3.4.1  Percent  Difference  Maps 

The  percent  difference  map  is  the  primary  method  of  data  representation  in  this 
paper.  To  create  the  difference,  two  TEC  outputs  from  GAIM-GM  must  be  generated.  A 
percent  difference  is  calculated  for  each  grid  point  in  the  GAIM-GM  output  using 
Equation  4. 
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TEC  output  A  and  B  represent  the  experiment  and  control  cases,  respectively.  For 
example,  to  show  the  effect  of  bubbles  on  the  TEC,  the  model  output  with  bubbles 
present  would  be  assigned  to  A  and  the  output  with  bubbles  absent  to  B.  Since  bubbles 
are  depleted  TEC  values  within  the  ionosphere,  a  negative  percent  difference  indicates 
the  presence  of  plasma  bubbles  in  the  GAIM-GM  output  because  B,  the  control  case  with 
bubbles  absent,  would  be  larger  than  A,  the  experiment  with  bubbles  present.  The 
percent  difference  output  is  then  plotted  against  the  latitude  and  longitude  axes  with  a 
world  map  overlay.  Percent  difference  is  a  more  meaningful  calculation  than  a  regular 
difference  between  the  TEC  outputs  because  the  percent  difference  contains  information 
about  the  control  state  of  the  ionosphere  and  how  much  deviation  the  experiment  state  of 
the  model  has  changed.  For  a  depletion  factor  of  10  the  percent  difference  from  the 
control  TEC  output  would  ideally  be  90%,  a  factor  of  100  would  be  99%,  and  a  factor  of 
1000  would  be  99.9%.  Percent  difference  maps  allow  for  a  quick  analysis  of  the 
accuracy  of  depletion  magnitude. 

3.4.2  Skill  Scores 

Skill  scores  were  utilized  as  a  method  to  objectively  evaluate  the  performance  of 
GAIM-GM.  A  skill  score  is  composed  of  three  inputs,  a  reference,  an  observation,  and  a 
tested  parameter.  The  reference  is  the  climatology,  in  this  case  the  GAIM-GM  output 
with  bubbles  absent.  The  IFM  TEC  output  with  plasma  bubbles  imposed  was  used  as  the 
observation  or  truth.  The  tested  parameter  was  the  GAIM-GM  output  constructed  with 
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GPS  input  files  having  bubbles  present.  A  skill  score  is  the  ratio  of  the  root  mean  square 
errors  of  the  inputs  given  by  Equations  5,  6,  and  7. 
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The  IFM  TEC  output  grid  and  the  GAIM-GM  TEC  output  grid  are  different 
resolutions.  In  order  to  use  the  skill  score  equations  the  grids  of  the  three  inputs  have  to 
match  exactly.  The  solution  is  to  tailor  the  GAIM-GM  regional  output  area  to  match  the 
IFM  grid  locations  by  adjusting  the  selected  region  in  the  GAIM-GM  set  up  file.  The 
grid  used  for  the  regional  GAIM-GM  output  and  the  skill  scores  goes  from  33S  to  10. 5N 
and  240E  to  326.25E.  The  resolution  of  the  output  is  1.5  degrees  of  latitude  by  3.75 
degrees  of  longitude,  the  same  resolution  as  the  modified  IFM  grid.  Skill  scores  range 
from  an  upper  limit  of  1  down  through  0  and  have  no  negative  bound.  A  skill  score 
greater  than  zero  means  the  tested  parameter’s  output  was  an  improvement  to  the 
reference  compared  to  the  observation.  A  negative  skill  score  means  the  tested 
parameter’s  output  was  worse  than  the  reference  compared  to  the  observation.  Skill 
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scores  closest  to  1  are  desired.  A  comparison  between  two  tested  parameters  can  occur  if 
the  same  reference  and  observation  are  used  to  compute  each  skill  score. 


3.5  Test  Cases 

GAIM-GM  was  run  with  and  without  bubbles  present  through  a  range  of  model 
parameters,  geophysical  conditions,  and  bubble  characteristics,  presented  in  Table  2. 
Model  parameters  tested  included  the  resolution,  time  constant,  and  ground  station 
density.  Both  the  global  and  regional  model  resolutions  were  tested  and  due  to  the  better 
performance  of  the  regional  mode,  most  of  the  remaining  test  cases  were  performed  in 
the  regional  resolution.  The  default  time  constant  of  5  hours  and  the  minimum  of  1  hour 
were  compared  in  both  global  and  regional  resolutions.  The  AFWA  GPS  ground  stations 
list  and  an  ideal  grid  of  regionally  distributed  ground  stations  were  used  to  test  the  effect 
of  varying  the  ground  station  density.  The  global  distribution  of  AFWA  ground  stations 
is  shown  in  Figure  9  and  regionally  in  Figure  10.  The  ideal  grid  of  regionally  distributed 
ground  stations  is  presented  in  Figure  11. 


Table  2.  List  of  parameters  varied  during  the  assessment  of  plasma  bubble  effects  on  GAIM-GM. 


Model  Parameters 

Geophysical  Conditions 

Bubble  Characteristics 

Model  Resolution 

1 )  Global 

2)  Regional 

Seasons 

1)  Spring  Equinox 

2)  Winter  Solstice 

Factor  of  Depletion 

1)  10 

2)  100 

3)  1000 

Time  Constant 

1 )  5  hours  (default) 

2)  1  hour  (minimum) 

Solar  Cycle 

1)  Solar  minimum 

F10.7  =  70  SFUs 

2)  Solar  medium 

F10.7  =  113  SFUs 

3)  Solar  maximum 

F10.7  =  210  SFUs 

Bubble  Movement 

1 )  Co-rotating 

0  m/s 

2)  Super-rotating 

230  m/s  to  the  east 

3)  Sun  synchronous 

660  m/s  to  the  west 

Ground  Station  Density 

1)  AFWA  Stations 

2)  Ideal  Grid 
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Figure  9.  Global  distribution  of  59  AFWA  GPS  ground  stations. 


Figure  10.  Regional  distribution  of  AFWA  station  list  GPS  ground  stations. 
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Figure  11.  Regional  distribution  of  ideal  GPS  ground  station  list. 

Next,  seasonal  and  solar  cycle  changes  were  made  to  test  geophysical  constraints 
on  plasma  bubble  reproduction  in  GAIM-GM.  The  spring  equinox  and  winter  solstice  of 
2004  were  chosen  for  the  seasonal  comparison.  Three  solar  cycle  conditions  were 
studied,  solar  minimum,  medium,  and  maximum.  The  three  solar  states  were 
characterized  by  their  F10.7  values,  which  are  a  measure  of  the  solar  flux  at  10.7  cm, 
which  is  used  as  a  proxy  for  the  shorter  wavelength  flux  that  produces  photoionization  in 
the  ionosphere.  Solar  minimum  was  defined  as  an  FI 0.7  =  70  solar  flux  units  (SFUs), 
solar  medium  as  a  F10.7  =113  SFUs,  and  solar  maximum  as  a  F10.7  =  210  SFUs. 
Geomagnetic  quiet  conditions,  Kp  =  1,  were  used  for  the  duration  of  the  research. 

Finally,  depletion  factors  and  bubble  movement  constituted  the  plasma  bubble 
characteristics  that  were  tested.  Three  factors  of  depletion,  10,  100,  and  1000  were  used. 
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Additionally,  three  methods  of  bubble  movement  were  studied,  co-rotating,  super- 
rotating,  and  sun  synchronous  bubble  movement.  Super-rotating  bubbles  moved  to  the 
east  at  230  m/s  after  their  formation  at  2000  LT.  Although  bubbles  have  only  been 
observed  to  super-rotate  at  speeds  up  to  190  m/s,  230  m/s  was  chosen  as  a  realistic  upper 
bound  that  fit  well  into  the  modified  IFM  grid  resolution.  The  sun  synchronous  bubble 
movement  is  physically  unrealistic,  but  provides  a  test  of  the  perturbation  movement 
assumption  built  into  the  Kalman  filter  by  moving  a  single  bubble  to  the  west  at  660  m/s, 
or  15  degrees  per  hour.  Table  3  lists  the  18  GAIM-GM  runs  that  were  conducted  by 
varying  these  parameters. 
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Table  3.  Variable  matrix  for  GAIM-GM  model  runs  completed. 


Run 

Number 

Day 

Mode 

Bubble  # 
or  type 

Depletion 

Factor 

Time 

Constant 

#  of  Ground 

Stations 

F10.7 

Remarks 

1 

82 

global 

0 

n/a 

5 

59 

113 

global  control 

2 

82 

global 

4  corotate 

10 

5 

59 

113 

global  bubbles  present 

3 

82 

regional 

0 

n/a 

5 

59 

113 

regional  control 

4 

82 

regional 

4  corotate 

10 

5 

59 

113 

regional  bubbles  present 

5 

82 

global 

4  corotate 

10 

1 

59 

113 

global  time  constant  =  1 

6 

82 

regional 

4  corotate 

10 

1 

59 

113 

regional  time  constant  =  1 

7 

356 

regional 

0 

n/a 

5 

59 

113 

regional  solstice  control 

8 

356 

regional 

4  corotate 

10 

5 

59 

113 

regional  solstice  factor  10  depletion 

9 

356 

regional 

4  corotate 

100 

5 

59 

113 

regional  solstice  factor  100  depletion 

10 

356 

regional 

4  corotate 

1000 

5 

59 

113 

regional  solstice  factor  1000  depletion 

11 

82 

global 

1  sun 

synchronous 

10 

5 

59 

113 

global  sun  synchronous  bubble 

12 

82 

regional 

4  super-rotate 

10 

5 

59 

113 

regional  super-rotate  bubbles 

13 

82 

regional 

4  corotate 

10 

5 

408 

113 

regional  ideal  grid  bubbles  present 

14 

82 

regional 

0 

n/a 

5 

408 

113 

regional  ideal  grid  control 

15 

82 

regional 

0 

n/a 

5 

59 

70 

solar  minimum  control 

16 

82 

regional 

4  corotate 

10 

5 

59 

70 

solar  minimum  bubbles  present 

17 

82 

regional 

0 

n/a 

5 

59 

210 

solar  maximum  control 

18 

82 

regional 

4  corotate 

10 

5 

59 

210 

solar  maximum  bubbles  present 

IV.  Results  and  Analysis 


This  chapter  will  present  and  analyze  the  results  of  placing  plasma  bubbles  into 
the  GPS  input  files  of  GAIM-GM.  First,  GAIM-GM  TEC  output  from  both  global  and 
regional  modes  will  be  compared  to  show  the  utility  of  the  regional  mode.  Next,  plasma 
bubbles  with  different  depletion  factors  will  be  evaluated  with  a  focus  on  the  similarities 
in  the  results.  Seasonal  and  solar  cycle  differences  will  follow  emphasizing  the  impact  of 
the  background  ionospheric  TEC  values  on  GAIM-GM  plasma  bubble  reconstruction. 
Then,  the  time  constant  model  parameter  will  be  evaluated  by  varying  the  time  constant 
and  bubble  motion.  Finally,  the  GPS  ground  station  density  will  be  increased  to  an  ideal 
regional  grid  to  showcase  the  effect  of  additional  ground  stations  on  plasma  bubble 
recreation. 

The  results  presented  in  this  thesis  expand  upon  an  earlier  study  on  the  effect  of 
plasma  bubbles  on  the  output  of  GAIM-GM  conducted  at  USU  by  the  Center  for 
Atmospheric  and  Space  Sciences.  Drs.  Schunk,  Scherliess,  and  Thompson  tested  the 
capability  of  GAIM-GM  to  model  plasma  bubbles  and  storm  enhanced  density  regions  in 
terms  of  the  quantity  of  GPS  observations  available  to  the  model  (2010).  Plasma  bubble 
characteristics  and  model  conditions  were  duplicated  from  the  USU  study  to  initialize  and 
verify  the  methods  and  models  used  presently.  Both  the  USU  and  present  methodology 
place  four  bubbles,  each  forming  at  2000  LT,  over  South  America  with  factor  of  10 
density  depletions.  The  only  significant  difference  between  the  two  procedures  was  the 
width  of  the  bubbles,  which  was  set  at  5  degrees  by  USU  and  7.5  in  this  paper  due  to  IFM 
grid  spacing.  The  separation  between  bubbles  in  both  cases  approximated  the  bubble 
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width  and  the  same  bubble  center  points  were  used.  Initial  output  for  this  study  was 
verified  against  the  USU  study  for  consistency  and  showed  the  similar  input  conditions 
independently  produce  comparable  results  in  GAIM-GM  output  (Schunk  et  al,  2010). 
The  verification  process  also  included  tests  of  the  time  constant,  latitude  correlation, 
plasmasphere  correction,  and  satellite  and  station  bias  correction,  all  of  which  performed 
as  expected. 

4.1  Global  and  Regional  Resolution  Analysis 

The  longitudinal  resolution  of  GAIM-GM  in  global  mode  is  15  degrees  which  is 
twice  as  large  as  the  7.5  degree  plasma  bubbles  used  in  this  study.  Because  the  model 
grid  resolution  is  larger  than  the  bubble  itself  it  was  anticipated  that  GAIM-GM  would 
have  difficulties  reproducing  the  bubbles,  particularly  the  magnitude  of  the  depletions. 
GAIM-GM  in  regional  mode  has  a  longitudinal  resolution  of  3.75  degrees,  half  that  of 
one  bubble  width,  and  therefore  was  expected  to  reproduce  bubbles  better  than  in  global 
mode.  Figure  12  shows  GAIM-GM  TEC  output  and  percent  difference  maps  for  both 
global  and  regional  model  runs.  Panels  (c)  and  (d)  are  GAIM-GM  TEC  output 
constructed  from  input  files  containing  plasma  bubbles  with  a  depletion  factor  of  10.  The 
four  plasma  bubble  structures  that  were  placed  into  the  input  files  cannot  be  resolved  in 
either  the  global  or  regional  resolution  TEC  outputs.  Panels  (e)  and  (f)  show  that  the 
depletions  present  in  the  TEC  outputs  do  not  resemble  the  bubble  positions  input. 
Additionally,  the  magnitude  of  the  depletions  is  less  than  30%  which  is  significantly  less 
than  the  90%  depletions  placed  in  the  input  files.  The  area  of  the  depletion  in  the  global 
mode  is  incorrectly  distributed  zonally,  seen  in  panel  (e).  The  bubbles,  as  input  into 
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Figure  12.  Comparison  of  GAIM-GM  global  and  regional  resolutions.  Panels  (a)  through  (d)  in 
units  of  TEC  *  1016  electrons/m2.  Panels  (e)  and  (f)  are  percent  differences  between  the  above  panels. 
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GAIM-GM  through  the  GPS  measurements,  are  elongated  north  to  south,  previously 
shown  in  Figure  6.  The  regional  mode  also  does  not  correctly  reproduce  the  bubbles,  but 
is  more  adept  at  keeping  the  bubbles’  longitudinal  distribution  than  the  global  mode.  The 
regional  mode  outperforms  the  global  mode  due  to  the  improved  ratio  of  grid  resolution 
to  feature  size.  Next,  regional  mode  was  used  to  study  the  effect  of  depletion  magnitude 
on  GAIM-GM  output. 

4.2  Magnitude  of  Depletion  Analysis 

A  comparison  of  depletion  factors  of  10,  100,  and  1000  is  presented  in  Figure  13, 
again  for  the  case  of  four  plasma  bubbles  over  South  America  at  03UT.  Panels  (a),  (b), 
and  (c)  all  look  remarkably  similar  and  in  fact  are  within  1.5  TECU,  or  10%.  The 
depletions  cover  a  nearly  identical  area  in  each  case.  Panels  (d)  and  (e)  confirm  that  all 
three  outputs  for  varying  depletion  factors  are  within  10%  of  each  other. 

The  observed  result  that  varying  the  magnitude  of  the  bubble  depletions  has  little 
affect  on  the  GAIM-GM  output  can  be  explained  by  an  analysis  of  the  Kalman  filter 
present  in  the  model.  The  Kalman  filter  in  GAIM-GM  discards  GPS  measurements 
greater  than  a  factor  of  5  less  than  the  background  ionosphere  created  from  the  previous 
time  step.  Because  different  look  angles  between  GPS  satellites  and  ground  stations  will 
see  varying  degrees  of  the  bubble,  not  all  look  angles  will  be  passed  through  the  filter. 
Figure  14  is  a  layout  of  three  example  look  angles  through  a  plasma  bubble  placed  in  the 
IFM  electron  density  grid.  Look  angle  1  does  not  pass  through  the  plasma  bubble  and  the 
associated  TEC  value  is  unaffected.  The  second  look  angle  passes  partially  through  the 
plasma  bubble.  TEC  values  resulted  from  look  angle  2  are  reduced,  but  not  by  enough  to 
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Figure  13.  Percent  difference  comparison  of  factor  of  10, 100,  and  1000  bubble  depletions.  Valid  at 
03UT. 
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Figure  14.  Selected  look  angles  through  a  plasma  bubble  as  viewed  in  the  modified  IFM  electron 
density  grid.  Look  angle  1  does  not  pass  through  the  bubble,  look  angle  2  partially  passes  through 
the  bubble,  and  look  angle  3  passes  through  the  entire  bubble. 


be  rejected  by  the  Kalman  filter.  Look  angle  2  produces  the  measurements  of  the  bubble 
that  are  assimilated  by  the  model  and  affected  by  the  depletion  magnitude.  Look  angle  3 
travels  through  most  of  the  bubble  and  its  TEC  values  are  reduced  beyond  a  depletion 
factor  of  5  resulting  in  a  rejection  of  the  measurement  by  the  Kalman  filter  for  all  three 
depletion  factors. 

Samples  of  elevation  angles  and  slant  TECs  are  given  in  Table  2  for  conditions 
where  there  were  depletion  factors  of  10,  100,  and  1000,  as  well  as  the  bubble  absent 
case.  The  three  elevation  angles  are  representative  of  the  look  angles  in  Figure  14.  For 
the  look  angle  1  case  of  17  degrees  the  slant  TEC  value  is  unaffected  by  the  depletion 
magnitude  because  that  look  angle  does  not  pass  through  the  bubble.  In  fact,  for  a  ground 
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Table  4.  Slant  TEC  values  measured  from  3  selected  look  angles  through  a  plasma  bubble  placed  in 
the  modified  IFM  electron  density  grid. 


No  bubbles 

lOx  depletion 

lOOx  depletion 

lOOOx  depletion 

Elevation  Angle 

Look  Angle  # 

Slant  TEC 

Slant  TEC 

Slant  TEC 

Slant  TEC 

17.23 

1 

46.93 

46.93 

46.93 

46.93 

43.02 

2 

30.66 

23.13 

22.48 

22.41 

63.55 

3 

21.15 

3.70 

1.97 

1.79 

station  centered  below  a  plasma  bubble  based  at  200  km  with  a  width  of  830  km  or  7.5 
degrees,  any  elevation  angle  less  than  25  degrees  will  not  pass  through  the  bubble.  Next, 
the  look  angle  2  case  of  43  degrees  passes  through  part  of  the  undisturbed  ionosphere  and 
part  of  the  bubble  and  is  therefore  only  partially  affected  by  the  bubble’s  presence.  The 
slant  TECs  for  depletion  factors  of  10,  100,  and  1000  are  reduced  by  factors  of  1.33,  1.36, 
and  1.37,  respectively.  None  of  the  depletion  factors  for  look  angle  2  are  greater  than  the 
factor  of  5  cutoff  of  the  Kalman  filter  and  therefore  the  measurement  will  be  used  in  all 
three  instances.  Although  the  factor  of  1000  depletion  produces  a  lower  slant  TEC  for 
the  elevation  angle  than  the  factor  of  10  depletion,  it  is  only  2%  lower  and  consequently 
cannot  have  a  large  affect  on  GAIM-GM’s  output.  The  effect  of  this  small  difference  in 
look  angle  2  cases  can  be  seen  in  panels  (d)  and  (e)  of  Figure  13.  Lastly,  the  high 
elevation  angle  of  63  degrees  has  depletion  factors  of  5.71,  10.74,  and  1 1.82,  all  of  which 
are  above  the  factor  of  5  threshold  of  the  Kalman  filter.  None  of  the  observations  for  the 
high  elevation  angle  are  assimilated  by  the  filter  and  therefore  have  no  affect  on  the 
model  output  for  any  depletion  magnitude. 
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Small  differences  that  are  present  in  panels  (d)  and  (e)  of  Figure  13  can  be  found 
in  regions  within  5  degrees  of  ground  station  locations.  The  differences  are  not  found 
directly  over  the  ground  stations  because  these  measurements  would  have  had  large 
elevation  angles  like  look  angle  case  3  and  been  removed  by  the  Kalman  filter  for  all 
three  depletion  cases.  Instead,  the  largest  differences  are  found  away  from  the  ground 
station  location  due  to  look  angle  2  cases  and  corresponding  ionospheric  pierce  points 
that  are  removed  from  the  ground  station.  Since  the  magnitude  of  the  plasma  bubble 
depletion  has  little  effect  on  the  output  of  GAIM-GM,  a  factor  of  10  depletion  was  used 
for  the  remaining  test  cases.  Next,  geophysical  conditions  will  be  varied  to  study  the 
effect  on  GAIM-GM  output  with  bubbles  present. 

4.3  Geophysical  Conditions  Analysis 

The  capability  of  GAIM-GM  to  model  plasma  bubbles  is  affected  by  the  TEC 
values  of  the  background  ionosphere.  GAIM-GM  produces  weaker  perturbation  fields  in 
the  presence  of  lower  background  TEC  values.  Low  TEC  values  of  1-2  TECU  are 
commonly  found  during  solar  minimum  nights,  particularly  during  the  solstices, 
diminishing  the  capability  of  GAIM-GM  to  add  perturbations  to  the  background  state. 
Two  geophysical  conditions  of  relatively  higher  and  lower  background  TEC  values  were 
evaluated,  resulting  in  seasonal  and  solar  cycle  comparisons. 

4.3.1  Equinox  and  Solstice  Comparison 

The  results  of  the  spring  equinox  and  winter  solstice  conditions  for  00,  03,  and  06 
UT  are  presented  in  Figure  15.  Both  scenarios  included  four  plasma  bubbles  with  factor 
of  10  density  depletions  placed  in  sequence  over  South  America.  Although  the  FI 0.7 
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Figure  15.  Comparison  of  percent  difference  maps  from  equinox  and  solstice. 
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value  was  1 13  SFUs  for  both  cases,  the  background  TEC  values  of  the  equinox  run  were 
higher  by  5-10  TECUs  than  the  solstice  run  because  of  different  solar  zenith  angles  and 
seasonal  differences  in  the  background  neutral  densities.  The  result  of  this  difference  in 
background  TEC  values  was  that  the  GAIM-GM  output  during  equinox  showed  stronger 
depletions  than  during  the  solstice.  In  general,  GAIM-GM  output  generates  a  10% 
greater  depletion  with  the  higher  background  levels  during  equinox.  The  regions  covered 
by  the  depletions  are  similar  because  identical  ground  station  to  satellite  look  angles  were 
used  for  both  cases.  With  bubble  absent  TEC  values  over  South  America  differing  by  as 
much  as  25%  from  equinox  to  solstice,  GAIM-GM  produces  up  to  a  10%  stronger 
depletion  in  the  region.  Peak  skill  scores  for  the  equinox  case  are  0.39  while  for  the 
solstice  case  the  best  scores  are  0.34  for  the  same  time  period.  The  skill  scores  verify  that 
GAIM-GM  has  a  better  performance  while  ingesting  plasma  bubble  observations  when 
the  background  TEC  values  are  higher. 

4.3.2  Solar  Minimum  and  Maximum  Comparison 

The  next  test  was  to  run  GAIM-GM  with  plasma  bubbles  present  during  both 
solar  minimum  and  solar  maximum  conditions.  Snap  shots  of  the  output  taken  at  00,  03, 
and  06  UT  are  shown  in  Figure  16.  The  solar  maximum  graphs  have  larger  depletions 
than  the  solar  minimum  case  for  all  time  steps  and  the  percent  differences  are  generally 
twice  as  large.  At  the  first  time  step,  panels  (a)  and  (b),  the  largest  percent  difference  for 
solar  minimum  is  -10%  while  for  solar  minimum  is  it  -25%.  The  perturbation  fields 
move  sun  synchronously,  become  superimposed,  and  sum  west  of  the  bubble  locations. 
By  the  06UT  panels,  the  largest  percent  differences  for  solar  minimum  and  maximum  are 
-36%  and  -54%,  respectively.  The  solar  maximum  case  has  larger  depletions  in  all  areas 
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Figure  16.  Percent  difference  comparison  between  TEC  output  for  solar  minimum  and  maximum. 
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of  the  South  American  region  because  the  background  TEC  values  in  the  nightside 
ionosphere  are  larger  than  during  solar  minimum.  The  background  TEC  values  during 
solar  maximum  can  be  up  to  an  order  of  magnitude  larger  in  the  20  to  24LT  period  near 
the  equator.  Since  the  area  of  the  South  American  region  covered  by  the  depletions  is  the 
same  for  both  cases  the  only  difference  between  the  GAIM-GM  outputs  is  the  percent 
decrease  at  each  grid  location.  The  solar  maximum  case  was  closer  to  the  ideal  depletion 
percentage  of  90%  and  therefore  is  a  better  solution.  The  skill  scores  verify  that  GAIM- 
GM  performed  better  during  the  solar  maximum  case,  with  a  peak  skill  score  of  0.29 
compared  to  0. 18  for  the  solar  minimum  case. 

4.4  Relaxation  Time  Analysis 

The  relaxation  time,  or  time  constant,  controls  the  length  of  time  GAIM-GM  will 
take  to  relax  a  perturbation  back  to  zero.  GAIM-GM  uses  a  sun  synchronous  movement 
of  the  perturbation  field  which  causes  differences  to  translate  west  while  relaxing  the 
density  field  to  the  background  state  over  time.  On  the  percent  difference  maps  this 
feature  manifests  itself  as  a  field  of  differences  propagating  west. 

4.4.1  Changes  to  the  Time  Constant 

Lowering  the  time  constant  forces  perturbations  to  relax  to  the  background  state 
over  a  shorter  time  period.  Because  the  perturbations  are  propagated  sun  synchronously 
by  the  model,  the  result  of  reducing  the  time  constant  from  5  hours  to  1  hour  is  to  reduce 
the  westward  extent  of  the  TEC  differences,  shown  in  Figure  17.  Sun  synchronous 
movement  means  that  a  perturbation  would  translate  15  degrees  per  hour  to  the  west. 
Reducing  the  time  constant  from  5  hours  to  1  hour  should  reduced  the  westward  extent  of 


47 


(a)  Time  constant  =  5  hours  (b)  Time  constant  =  1  hour 

Figure  17.  Comparison  of  time  constant  settings  for  08UT  on  day  82. 


differences  by  60  degrees.  Figure  17  shows  a  reduction  of  the  westward  extent  by  45 
degrees  and  in  general  reducing  the  time  constant  reduced  the  westward  limit  of  the 
depletions  between  30  and  60  degrees. 

In  addition  to  changing  the  westward  boundary  of  the  perturbation  field,  Figure 
17  also  shows  that  the  depletion  had  a  smaller  magnitude  when  the  time  constant  was  set 
to  one.  Lowering  the  time  constant  not  only  affects  the  distance  to  the  west  a  disturbance 
can  propagate,  but  also  the  rate  at  which  it  decays  to  the  background  state.  A  shorter 
relaxation  time  means  the  perturbation  has  to  decay  more  quickly  and  by  a  larger 
magnitude  at  each  time  step.  Four  bubbles  were  used  in  this  test,  and  each  created  their 
own  perturbation  fields  that  summed  as  the  fields  became  superimposed  upon  one 
another.  The  result  of  the  summed  perturbation  fields  was  to  create  the  largest  magnitude 
depletions  west  of  South  America.  The  lower  time  constant  meant  each  of  the  four 
perturbations  decayed  faster  resulting  in  a  smaller  magnitude  depletion  than  in  the  default 
time  constant  case. 
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Another  factor  affecting  the  perturbation  fields  is  the  distribution  of  ground 
stations.  When  the  perturbation  is  propagated  into  a  region  that  has  GPS  measurements 
for  that  time  step,  the  presence  of  measurements  from  the  ground  station  causes  the 
perturbation  field  to  quickly  relax.  The  sun  synchronous  extent  of  any  perturbation  field 
is  therefore  shaped  by  the  distribution  of  ground  stations  to  the  west  of  the  initial 
perturbation. 

While  a  smaller  spatial  extent  is  affected  by  the  perturbation  field  when  a  lower 
time  constant  is  used,  the  skills  scores  are  actually  better  for  the  higher  time  constant. 
Figure  18  shows  the  skill  scores  for  regional  runs  of  GAIM-GM  with  time  constants  set 
to  five  and  one  hour.  The  higher  skill  scores  with  the  time  constant  set  to  five  hours 
indicate  that  the  output  was  closer  to  the  ideal  observation  set,  the  IFM  with  bubbles 
present,  than  the  lower  time  constant  case.  The  reason  for  this  result  is  that  while  the 
smaller  disturbed  area,  resulting  from  the  time  constant  being  set  to  one  hour,  initially 
appears  better,  it  is  actually  the  larger  magnitude  depletions  present  off  the  west  coast  of 
South  America  in  the  time  constant  set  to  five  hour  case  that  come  closer  to  the  actual 
bubble  observations.  There  is  a  balance  between  the  spatial  extent  and  magnitude  of  the 
depletion  in  the  GAIM-GM  TEC  output.  Lowering  the  time  constant  improves  the 
spatial  boundaries  but  degrades  the  reproduction  of  depletion  magnitude,  yielding  a  net 
worse  result. 

4.4.2  Sun  Synchronous  Bubble  Movement 

Since  GAIM-GM  moves  perturbations  sun  synchronously  the  next  test  case  was 
to  allow  the  bubble  to  move  in  this  manner.  Although  this  condition  is  physically 
unrealistic,  it  provides  an  interesting  result  by  aligning  bubble  movement  with  the  sun 
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Figure  18.  Skill  score  comparison  of  varying  the  time  constant  in  GAIM-GM  regional  output  for  day 
82. 

synchronous  assumption  within  the  Kalman  filter.  One  bubble  was  created  over  Brazil 
and  allowed  to  move  sun  synchronously  15  degrees  per  hour  to  the  west.  The  result  of 
this  run  is  presented  in  Figure  19  as  a  series  of  percent  difference  graphs  between  the  sun 
synchronous  case  and  a  run  of  GAIM-GM  with  plasma  bubbles  absent.  Panels  (a) 
through  (e)  are  snapshots  of  the  GAIM-GM  output  at  one  hour  increments  from  00-04 
UT.  The  bubble  first  appears  at  00  UT  and  the  area  and  magnitude  of  the  depletion  in  the 
GAIM-GM  output  is  similar  to  the  00  UT  frame  from  previous  comparisons.  As  time 
progresses  in  panels  (b)  through  (e)  the  perturbation  field  moves  to  the  west  at  the  same 
rate  as  the  bubble.  The  measurements  GAIM-GM  is  receiving  through  the  GPS  inputs 
are  matching  the  location  of  the  perturbation  field.  The  perturbation  field  reinforces  the 
measurements  at  each  time  step  resulting  in  the  model  successfully  tracking  the  location 
of  the  bubble  through  04  UT.  The  width  of  the  bubble  is  approximately  1 1 .25  degrees, 
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Figure  19.  Percent  difference  maps  for  GAIM-GM  regional  output  on  day  82  with  sun  synchronous 
bubble  measurements  as  inputs  to  the  model. 
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about  one  and  a  half  times  larger  than  the  7.5  degree  bubble  that  the  observations  were 
created  from.  The  depletion  magnitude  is  underrepresented  by  GAIM-GM,  with  a 
maximum  depletion  of  35%  being  significantly  less  than  the  90%  depletion  that  was 
input  into  the  GPS  measurements. 

4.4.3  Super- Rotating  Bubble  Movement 

Plasma  bubble  super-rotation  has  been  observed  and  was  tested  by  moving  the 
bubbles  to  the  east  at  230  m/s,  or  7.5  degrees  per  hour,  after  their  formation  at  20LT.  The 
results  are  displayed  in  Figure  20  as  percent  difference  maps.  As  expected  from  the  sun 
synchronous  assumption  in  the  Kalman  filter,  the  model  was  unable  to  move  the 
perturbation  fields  in  the  correct  direction  towards  the  east  and  instead  moved  the 
depletions  sun  synchronously  to  the  west.  The  perturbations  fields  were  relaxed  when 
bubble  absent  measurements  over  the  Galapagos  and  Easter  Islands  were  encountered, 
but  extended  farther  west  in  regions  that  ground  stations  were  not  present.  At  no  point 
was  the  model  able  to  move  the  perturbation  field  into  the  Atlantic  east  of  the  ground 
station  locations.  Figure  7  showed  the  movement  of  the  bubbles  in  the  IFM  grid. 

4.5  Station  Density  Analysis 

GAIM-GM  output  with  the  AFWA  station  list  is  compared  to  output  with  an 
ideal,  regional  grid  of  GPS  ground  stations  in  Figure  21.  The  ideal  grid  of  ground 
stations  provides  a  dense  network  of  GPS  measurements  to  the  Kalman  filter  from  which 
GAIM-GM  is  able  to  reproduce  four  distinct  bubbles.  The  density  depletions  are 
centered  at  the  correct  longitudes  to  match  the  bubble  input  placement  but  the  depletion 
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Figure  20.  Percent  difference  maps  for  GAIM-GM  regional  output  on  day  82  with  super-rotating 
bubble  measurements  as  inputs  to  the  model. 
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Figure  21.  Percent  difference  map  of  4  corotating  plasma  bubbles  in  an  ideal  regional  grid  for  day 
82. 

magnitudes  are  below  the  bubble  measurements.  Depletion  magnitudes  in  the  GAIM- 
GM  output  were  35-40%,  significantly  less  than  the  90%  depletion  used  in  the  input  files. 

The  output  with  the  ideal  grid  supports  results  from  the  time  constant  and 
geophysical  conditions  tests.  First,  GAIM-GM  is  able  to  create  individual  bubbles 
because  the  ideal  grid  density  is  great  enough  to  place  a  line  of  ground  stations  between 
each  bubble.  The  perturbation  field  created  from  the  measurements  at  each  bubble  is 
propagated  sun  synchronously,  but  as  discussed  above  is  translated  to  a  point  with 
measurements  from  a  new  ground  station.  The  measurements  with  bubbles  absent  from 
the  ground  station  to  the  west  are  weighted  more  heavily  in  the  Kalman  filter  and  the 
perturbation  field  is  relaxed  to  zero.  Second,  the  ideal  grid  output  has  lower  depletion 
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percent  differences  in  the  eastern  bubble  as  compared  to  the  western  bubble.  The  eastern 
bubble  lies  in  an  area  of  lower  background  TEC  values  than  the  western  bubble  and 
supports  the  earlier  data  that  GAIM-GM  more  reliably  recreates  bubble  observations  in 
higher  TEC  environments. 
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V.  Discussion 


The  results  of  the  plasma  bubble  output  analyses  are  discussed  in  this  chapter, 
including  their  computational  and  operational  uses  and  applications.  Additionally,  topics 
for  further  research  will  be  suggested  and  a  brief  summary  will  follow. 

5.1  GAIM-GM  Output  Analysis 

GAIM-GM  is  limited  in  its  ability  to  model  plasma  bubbles  because  of  several 
parameters  within  the  Kalman  filter.  First,  the  lower  bound  the  filter  places  on  acceptable 
observed  TEC  values  means  that  the  highest  magnitude  depletion  measurements  will  be 
rejected  and  not  used  within  the  GAIM-GM  calculation.  GPS  measurements  below  this 
threshold  are  discarded  by  the  filter  because  they  are  too  far  below  the  background  model 
state.  While  it  is  important  to  protect  against  erroneous  data,  in  this  case  the  outlier  data 
represents  a  real  physical  feature  of  the  ionosphere.  Because  the  filter  removes 
observations  below  a  factor  of  5  depletion  from  the  background  state,  GAIM-GM  cannot 
depict  the  full  magnitude  of  plasma  bubbles  with  depletion  factors  of  5  and  greater. 

The  next  limiting  factor  within  the  Kalman  filter  is  the  assumption  that 
perturbations  move  sun  synchronously.  Since  plasma  bubbles  have  been  observed  to 
corotate  with  the  Earth  and  even  super-rotate  this  assumption  is  not  valid  for  this 
particular  kind  of  perturbation.  The  result  of  the  sun  synchronous  assumption  is  to 
incorrectly  move  the  perturbation  field  to  the  west  at  a  rate  of  15  degrees  per  hour. 
Because  the  default  time  constant  within  the  model  is  5  hours  this  can  result  in  the 
perturbation  field  extending  more  than  75  degrees  to  the  west  in  the  absence  of  GPS 
measurements.  The  perturbation  field  is  erroneously  changing  the  output  TEC  values 
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west  of  the  actual  plasma  bubble  location.  In  this  study  TEC  changes  were  documented 
up  to  100  degrees  to  the  west  for  a  period  of  up  to  5  hours  after  the  plasma  bubble.  This 
study  evaluated  the  prospect  of  lowering  the  time  constant,  which  was  shown  to  limit  the 
sun  synchronous  movement  of  the  perturbation  field,  however  the  effect  of  this  change  on 
other  features  of  the  ionosphere  was  not  accounted  for. 

Operationally,  GAIM-GM’s  inability  to  capture  the  full  magnitude  of  bubble 
depletions  and  the  false  depletion  area  to  the  west  of  the  bubble  location  can  lead  to 
incorrect  specifications  and  forecasts  of  the  ionosphere.  The  electron  density  constructed 
by  GAIM-GM  will  be  incorrect  in  these  areas  leading  to  errors  in  the  TEC  and  height  and 
density  of  the  F2  layer  peak.  The  incorrect  specifications  and  forecasts  will  lead  to 
incorrect  range  and  frequency  specifications  for  HF  radio  communication. 

Coupled  to  any  perturbation  field  discussion  is  the  dependence  of  GAIM-GM  on 
the  number  of  GPS  ground  stations  available.  In  the  absence  of  new  observations  the 
perturbation  field  is  imposed  upon  the  background  state.  In  the  presence  of  new 
observations  the  perturbation  field  is  changed  according  to  the  relative  weights  of  the 
new  observations.  This  result  explains  the  station  density  analysis  where  distinct  bubbles 
were  modeled  in  the  ideal  regional  grid.  Although  the  AFWA  station  list  of  59  stations 
may  seem  small  there  aren’t  that  many  more  equatorial  ground  stations  currently  in 
operation  that  could  be  added.  The  equatorial  regions  of  the  Earth  are  mostly  ocean  and 
low  population  density  land,  limiting  the  opportunity  for  a  dense  network  of  GPS  ground 
stations  such  as  is  the  case  of  the  continental  United  States  and  Europe.  Incorporating 
additional  satellite  data,  observing  the  ionosphere  near  the  dusk  sector,  and  using  GPS 
data  from  ships,  can  help  to  circumvent  some  of  these  problems. 
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The  resolution  of  GAIM-GM  is  also  a  limiting  factor  because  of  the  small  scale  of 
plasma  bubbles.  The  7.5  degree  longitudinal  width  of  the  bubbles  used  in  this  study  is 
reasonable,  but  extremely  large  for  a  bubble  and  was  used  because  it  fit  well  with  IFM 
resolution  while  still  being  realistic.  The  15  degree  longitudinal  global  resolution  of 
GAIM-GM  means  that  individual  bubbles  cannot  be  resolved.  Additionally,  any 
observations  of  plasma  bubble  depletions  will  be  averaged  and  smoothed  with  data  from 
a  large  region,  further  reducing  the  magnitude  of  the  depletion  that  GAIM-GM  can 
model. 

GAIM-GM  TEC  output  does  not  provide  sufficient  contrast  for  the  forecaster  to 
locate  plasma  bubbles.  Factor  of  10  bubble  depletions  are  only  1-2  TECU  differences 
from  the  background  state  in  the  nightside  ionosphere  which  would  be  difficult  to  find 
with  a  perfect  replication  of  the  bubble.  The  small  difference  coupled  with  the  diffuse 
bubble  output  from  GAIM-GM  makes  it  impractical  to  find  bubbles  in  the  TEC  output 
from  GAIM-GM.  Rather,  the  only  way  to  know  GAIM-GM  has  added  a  perturbation 
field  is  to  use  a  difference  map  such  as  the  ones  in  this  study.  This  would  require  two 
simultaneous  runs  of  GAIM-GM  to  compare  and  a  priori  knowledge  of  which  GPS  data 
files  contained  bubble  observations.  Even  after  creating  difference  maps  the  result  would 
not  be  much  more  useful  because  of  the  limitations  involving  the  Kalman  filter  discussed 
above. 

5.2  Further  Study 

Next,  options  for  further  study  are  suggested  that  would  further  test  the 
capabilities  of  GAIM-GM  to  model  plasma  bubbles. 
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5.2.1  Additional  Data  Sources 


The  present  study  only  incorporated  one  of  the  several  types  of  data  GAIM-GM  is 
capable  of  ingesting.  Satellite  occultations  could  be  used  as  an  additional  data  source  for 
the  model  rather  than  exclusively  GPS  slant  TEC  measurements.  The  advantage  of 
satellite  occultations  is  that  they  can  occur  in  any  geographic  location  and  are  not  tied  to 
ground  station  placement  like  GPS  measurements.  Occultation  measurements  would 
have  a  more  uniform  distribution  than  GPS  ground  stations  and  include  readings  over  the 
oceans  where  GPS  is  lacking.  Additionally,  incorporating  GPS  data  from  ships  would 
increase  the  density  of  GPS  measurements.  Many  routine  shipping  lanes  cross  equatorial 
regions  and  could  potentially  provide  a  wealth  of  additional  information.  Since  real  time 
measurements  outweigh  perturbation  fields  within  the  Kalman  filter  the  presence  of  a 
more  uniform  grid  and  data  over  the  oceans  should  limit  the  sun  synchronous  movement 
of  plasma  bubble  perturbation  field  yielding  a  better  result  than  GPS  ground  stations 
alone. 

5.2.2  Effect  of  Kalman  Filter  Modifications 

One  finding  of  the  present  paper  was  that  lowering  the  time  constant  produced  a 
spatially  improved  GAIM-GM  output  with  plasma  bubbles  because  the  perturbation 
fields  were  relaxed  to  background  values  faster.  Although  this  result  is  true  for  bubbles  it 
remains  unknown  how  this  change  affects  other  features  of  the  ionosphere.  Before  the 
time  constant  can  be  lowered  in  the  operational  model  the  effects  of  this  change  must  be 
investigated  to  understand  the  changes  in  other  regions  of  the  ionosphere  that  the  change 
will  invoke. 
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In  addition,  removing  the  sun  synchronous  movement  of  perturbation  fields 
should  result  in  a  better  plasma  bubble  output  from  GAIM-GM.  A  future  study  to 
investigate  this  alternative,  as  well  as  its  effects  on  other  features  of  the  ionosphere  could 
be  undertaken. 

5.3  Conclusion 

Overall,  GAIM-GM  is  not  capable  of  modeling  equatorial  plasma  bubbles.  The 
model  does  not  include  the  physics  necessary  to  account  for  bubbles  and  the  structure  of 
the  Kalman  filter  results  in  an  inaccurate  recreation  of  the  local  ionosphere.  Although  the 
model  is  able  to  lower  the  TEC  values  in  the  region  near  the  bubble,  the  area  of  that 
depletion  is  too  large  and  the  magnitude  of  the  depleted  values  is  too  small.  Moreover, 
the  sun  synchronous  movement  of  perturbations  by  the  Kalman  filter  propagates  depleted 
values  into  regions  where  bubbles  are  not  present. 

The  research  presented  in  this  paper  is  important  to  both  the  model  developers 
and  the  Air  Force.  By  understanding  the  limitations  of  the  data  assimilation  scheme, 
ionospheric  model  developers  will  be  able  to  create  improved  products  in  the  future. 
Additionally,  the  Air  Force  can  use  this  research  to  understand  how  the  presence  of 
plasma  bubbles  will  affect  the  output  of  its  primary  ionospheric  model.  Knowledge  of 
plasma  bubble  effects  on  GAIM-GM  output  will  enable  forecasters  to  improve  forecasts 
and  tailor  products  to  customer  needs.  Decision  makers  will  have  a  clearer  picture  of  the 
capabilities  and  limitations  of  GAIM-GM  leading  to  a  greater  probability  of  mission 
success. 
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